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mid-infrared spectral analysis methods 

A. Juliasz^, Th. Henning^'^, J. Bouwman^, CP. Dullemond^, I. Pascucci^, D. Apai^'^ 

j uhaszOmpia-hd . mpg . de 

ABSTRACT 

The spectral region around 10 /zm, showing prominent emission bands from 
various dust species is commonly used for the evaluation of the chemical com- 
position of protoplanetary dust. Different analysis methods have been proposed 
for this purpose, but so far, no comparative test has been performed to test the 
validity of their assumptions. In this paper, we evaluate how good the various 
methods are in deriving the chemical composition of dust grains from infrared 
spectroscopy. 

Synthetic spectra of disk models with different geometries and central sources 
were calculated, using a 2D radiative transfer code. These spectra were then fitted 
in a blind test by four spectral decomposition methods. We studied the effect of 
disk structure (flared vs. flat), inclination angle, size of an inner disk hole and 
stellar luminosity on the fitted chemical composition. 

Our results show that the dust parameters obtained by all methods deviate 
systematically from the input data of the synthetic spectra. The dust composi- 
tion fitted by the new two-layer temperature distribution method, described in 
this paper, differs the least from the input dust composition and the results show 
the weakest systematic effects. The reason for the deviations of the results given 
by the previously used methods lies in their simplifying assumptions. Due to the 
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radial extent of the 10 /xm emitting region there is dust at different temperatures 
contributing to the flux in the sihcate feature. Therefore, the assumption of a 
single averaged grain temperature can be a strong limitation of the previously 
used methods. The continuum below the feature can consist of multiple com- 
ponents (e.g., star, inner rim, disk midplane), which cannot simply be described 
by a Planck function at a single temperature. In addition, the optically thin 
emission of "featureless" grains (e.g., carbon in the considered wavelength range) 
produces a degeneracy in the models with the optically thick emission of the disk. 

The effect of different noise levels on the results has also been tested. We 
find that for a singal-to-noise ratio of 100 one can expect an absolute uncertainty 
in the value of the crystallinity of about 11% using ground based observations 
(8-13 yum). For space based observations (7-17 /im) the expected uncertainty is 
about 5 % for the same signal-to-noise value. Moreover, the average value of the 
estimated crystallinity increases toward lower signal-to-noise ratios in general. On 
the basis of our results, we propose a recipe for the analysis and interpretation of 
dust spectroscopy data in the mid-infrared which should be especially valuable 
for analyzing Spitzer spectroscopy data and ground-based infrared spectroscopy 
data in the 10 yum window. 

Subject headings: astrochemistry - infrared: general - stars: circumstellar matter 
- protoplanetary disks - lines:profiles 



1. Introduction 



Excess emission above the stellar photosphere at mid-infrared wavelengths is a character- 
istic feature of Young Stellar Objects (YSOs). This infrared excess is caused by the thermal 
emission of circumstellar dust around the central star gathered in a disk-like structure and/or 
envelope. The amount of excess emission and the shape of the Spectral Energy Distrib ution 



(SED) is indicative of the evolutionary state of the object (e.g.. Ivan Boekel et al.ll2006l ). Not 



only geometrical and density the structure of the circumstellar material, b ut also the dust 



prop e rties evolve with time , due to, e.g., crystallization and coagulation fe.g.. lBeckwith et al. 



2000 : iHenning et al.l l2006l and references therein), vaporization and re-condensation (IGail 
ioM)- 



The dust evolution and grain processing can be even more clearly seen in the emis- 
sion/absorption features observed in the mid-infrared spectrum of many young stars with a 
wide range of spectral types (e.g., 
20081 : Henning 2008 (in press)). 



Molster fc Watersll2003l : lApai et al.ll2005l : iBouwman et al. 



The most frequently observed features, around 10 and 
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18 yum, are caused by amorphous and crystalline silicate grains. In the case of disk emission 
spectra, they arise from optically thin regions. Therefore, their shape is directly determined 
by the optical properties of the individual dust grains. The silicate features have been 
used in order to determine the dust properties, e.g., the fraction of crystalline grains, the 



iron / magnesium ratio or the characteristic grain s i ze in sur : 



(e.g., 


Waelkens et al. 


iggd; 


Bouwman et al. 


2( 


)01; 


Li et al. 


2004; 


van Boeke 


et a . 


2005b) or 


of Solar System comets ( 


Crovisier et al.l 


1997; 


Brucato et al. 


1999; 


Wooden 


2002; 


Min et al. 



ace layers of protoplanetary di sks 



2005bl ). Amorphous silicates with olivine and pyroxene stoichiometry are c ommon materi 



als in the diffuse interstellar medium (ISM) as well as in cicumstellar disks (IHenning et al. 



20051 ). The broad and smooth features with a peak at approximately 9.7 fim and 18.5 /im 
are attributed to the Si-0 stretching and 0-Si-O bending modes, respectively, in amor- 
phous silicate grains. Crystalline silicates (e.g., forsterite and enstatite) have in particular 
been observed in pr otoplanetary disks and they seem to be absent from the diffuse ISM. 
Kemper et al.l (120051 ) placed an 2.2% upper limit for the crystallinity fraction of the sili- 
cates in the diffuse ISM. Thus it seems to be reasonable to conclude that crystalline silicates 
form in protoplanetary disks. However, the details of the formation process remain un clear. 
Cryst alline silicate grains can be produced by direct condensation f rom the gas phase ( 



2001 



2004 ) or by annealing of ana orphous grains at higher temperatures (iFabian et al.l 12000 



Gail 



Gail 



Harker &: Deschll2002l ). or both processes. 



Space-based observations of disks with the Infrared Space Observatory (ISO) and the 
Spitzer Space Telescope (SST) as well as ground-based measurements with, e.g., TIMMI2, 
TReCS or COMICS, enabled the detailed investigation of the silicate features. Such stud- 



with different spectral types (e.g.. 


Bouwman et al. 


2001; 


van Boekel et a 


2005b; 


Apai et al. 


2005 


; iKessler-Silacci et a . 


2006; ] 


ionda et al. 


2006; 


Sareent et al. 


2006; ^ 


kVatson et al. 


2007; 


Sicilia-Aguilar et al. 


2007; 


Bouwman et al. 


2008 


). These studies also demonstrated that evi- 



dence for grain growth can be found in the mid-infrared spectra, indicated by the broadening 
of the 10 /im emission feature together with a decrease of its height above the underlying 
continuum. The relation, however, between grain growth and crystallization and the role 
of the stellar type is a topic of ongoing debat e. There are reports on possible correla- 



tions between grain grow th and crystallization (Ivan Boekel et al.l l2005bl ; lApai et al.l 12005 



Kessler-Silacci et al.ll2006r) and on an i nverse correlation between the average grain size and 
stellar type (IKessler-Silacci et al.ll2006l) and a possible inverse correlation between the crys- 
tallinity and the stellar temperature (lApai et al.ll2005l ). 



The results of these studies depend not only on the quality of observational data, but 
also on the applied analytical methods. Here, we should also note, that the features are 
produced in the surface layer and the material properties are not necessary typical of the 
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material in the disk midplane. In addition, the region emitting the 10 /xm feature will be 
located at different distances from the central star, depending on the radiation characteristics 
of the star and the nature of the dust grains and the disk properties. Due to the lack of 
a comparative study of the various analysis methods, it is practically impossible to find 
out if the relations, discussed above, are really present or if they are caused by systematic 
uncertainties in the applied methods. Since we can determine the physical properties and 
chemical composition of protoplanetary dust grains for a larger sample only via such spectral 
decomposition methods, it is important to know how good these methods are and which of 
them should be applied in a specific case. In this paper, using a 2D radiative transfer model 
of protoplanetary disks with a prescribed dust composition, we synthesize "observations" 
which we fit with three widely used spectral decomposition methods. We then compare the 
resulting dust composition to the original input composition. We also propose a new and 
fast method for fitting the mid-infrared spectral features in which we apply a distribution of 
temperatures, leading to more reliable results. 

2. Modelling 
2.1. Disk models 



In order to investigate the quahty of the spectral decomposition methods, we fitted syn- 
theti c spectra calculated by the 2D radiative transfer code RADMC (iDuUemond fc Dominik 
20041 ). Although this code assumes an axisymmetric density distribution, photon packages 
are followed in three dimensions by a Monte Carlo method. After the temperature distri- 
bution has been determined , the ray tracing module of the more general code RADICAL 
( jDuUemond fc Turollall2000l ) has been used to compute the SEDs. These codes have already 
been used to model protoplanetary disk SEDs in numerous s tudies ( iPontoppidan fc DuUemond 
20051 : Ivan Boekel et al.ll2005al : DuUemond fc pominikll2004r) and tested against other contin- 
uum radiative transfer codes (iPascucci et al. 1 l2004h . 



The disk structure can possibly be related to grain evolution. A flared disk, con- 
taining unprocessed, small amorphous grains becomes flatten ed with time as dust grains 
grow and settle to the disk midplane (iBouwman et al.ll2008l ). Since the flaring angle of 
the disk has a strong effect on the shape of the continuum at mid-infrared wavelengths 
( iDullemond fc Dominikll2004l ). it can also have an effect on the results of the spectral anal- 
ysis methods. Therefore, we used disk models with different geometries. We want to note 
that we did not treat the dust coagulation and sedimentation in a self-consistent way during 
the modelling. Since the irradiated area of a circumstellar disk depends on the luminosity of 
the central star, we used stars with three different types (Herbig Ae star, T Tauri star and 
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Brown Dwarf) as central radiation source in our models. The main parameters of the 2D 
RT disk models as well as the applied stellar properties are summarized in Tab.[Tl The stel- 
l ar emission h as been modelled by a blackbody radiation source instead of a Kurucz model 



( iKuruczl I1993I ). which made the models more uniform enabling us to compare the results 



directly. 

We calculated six series of disk models in order to investigate the effect of disk structure 
and disk orientation on the results of the fitting procedure. Each of the model series contained 
12 spectra: three different stellar types and four different inclinations from pole-on to close 
to edge-on orientation (0°, 20°, 45°, 70°). The mass of the disk in each model has been 
scaled with the stellar mass as Mdisk = O.OSACt. The surface density as a function of radius 
was described by a power law as S(r) = So('^/''^disk)^5 where So is the surface density at the 
outer disk radius and p = —1. The outer radius of the Herbig Ae disk was set to 400 AU. 
Then the outer disk radius of the T Tauri and Brown Dwarf disks was determined by the 
radial optical depth at 0.55 /im, which we took to be same as in the Herbig Ae disk. It is 
important to note that the chemical composition of dust grains was assumed to be uniform 
throughout the disk. This simplifying assumption will be less justified when considering a 
large wavelength range and, therefor e, a wider tempe r ature and radial region in the disk. 



The passive irradiated disk model of iDuUemond et al.l (1200 ll ) has been applied for the disk 



structure in general, but the detailed structure of the disks were different in the different 
model series. The six model series can be summarized as follows. 

Model AS, AB. A flared disk geometry has been applied with a flaring index a 
(see Eg. I A3 P of 2/7 and with a dimensionless pressure scale height at the inner radius of 
H{Rin) / Rin = 0.03 (see AppendixlA|. The inner disk radius was determined by the subli- 
mation temperature of 1500 K (Model AS) and 1060 K (Model AB). The models with 1060 K 
temperature at their inner edge have twice as large inner holes as the models with 1500 K 
temperature. 

Model BS, BB. Here we applied a moderately flared disk with a flare index of 1/7 
and with a dimensionless pressure scale height at the inner radius of H{Rin)/Rin = 0.03. 
The inner disk radius is again determined by the grain temperature of 1500 K (Model BS) 
and 1060 K (Model BB). 

Model CS, CB. A more flattened disk model has been considered here compared to 
the Model B series. We used a similar density structure for these models as it was used 
for Model B disks. For Model C the flaring index was set to 1/8 producing a self-shadowed 
flattened disk structure. The inner disk radius was again set by the grain temperature of 
1500 K (Model CS) and 1060 K (Model CB). 
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2.2. Mass absorption coefficients 

During the tests, we used six different dust species (amorplious silicates witli olivine 
and pyroxene stoichiometry, forsterite, enstatite, silica and carbon), each with three grain 
sizes (0.1, 1.5 and 6.0 /xm). Tab. [2] summarizes the used dust species, grain model and origin 
of the optical constants. The relative mass fraction for each dust species and grain size 
used in the study is presented in Tab. [31 The used relative mass fractions as a function of 
grain size corresponds to number density distribution of n oc a~^'^ and n oc a~^'^ for the 
amorphous and crystalline dust species, respectively. The mass absorption coefficients of 
these dust species as a function of wavelength are shown in Fig.[Tl In order to calculate the 
mass absorption coeffici ent from the opti cal constants, we used the theory of Distribution of 



Hollow Spheres (DHS) (IMin et al.ll2005al ) for crystalline silicates, while classical Mie theory 



for spherical particles was used for the amorphous dust species. 

There are two distinct, but coupled problems related to the fitting of the mid-infrared 
silicate features. One of the problems is the apphcability of a special set of optical data and 
the other one is the uncertainty in the applied spectral decomposition methods, themselves. 
It is unknown whether or not we use exactly the right optical data, which determine how 
the protoplanetary dust grains interact with radiation. The question of how good our set 
of mass absorption coefficients is is beyond the scope of this paper. Here, we used the 
same dust composition for all disk models (see Tab.E]), in order to evaluate the performance 
of the various spectral decomposition methods. We are aware that the use of one specific 
dust composition in all disk models is a limitation of the study, but it is impractical to 
test the methods on all possible combinations of dust mixtures and stellar/disk parameters. 
We chose, therefore, a single dust composition for the disk models the calculated spectra 
of which resembles that of typically seen towards young stellar objects. In this study we 
investigate how the different stellar and disk parameters can modify the results of spectral 
decomposition. The best way to isolate these effects is to use one dust mixture in all disk 
models, thus any difference between the input and the resulting dust parameters is certainly 
caused by the effects of stellar or disk parameters, which is then a deficiency of the applied 
method. 



2.3. Fitting methods 



Numerous techn iques have been applied for the spectral decomposition of mid-infrared 
silica t e features (e.g. , 



Bouwman et al.ll2001 



20061 : IWatson et aPboO?! : 



Bouwman et al. 



Li et al.ll2004 ; van Boekel et al.ll2005b ; Kessler-Silacci et al 



20081 ). but no systematic comparison of various 



techniques has been performed so far, leaving uncertainties in the quality and applicability 
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of these techniques. We choose the three most widely used methods and tested them against 
a new method described in this paper and apphed them to the synthetic spectra of the above 
described disk models. These three methods have been applied by many authors to evaluate 



wide range of source types (e.g., 


Bou 


wman et al. 


2001; 


Meeus et al. 


2003; 


van Boeke 


et al. 


2005b: 


Apai et al. 


2005 


; I 


ionda et al. 


2006 


Scheeerer et al. 


20od; 


Kessler-Silacci et al. 


2006; 


Sicilia-Aguilar et al. 


2007 


: Bouwman et al. 


2008). 


All methods assume that the silicate emis- 



sion comes from an optically thin region, but the temperature structure and underlying 
continuum is described in different ways in the distinct methods. 

The simplest method for the spectral decomposition of the silicate emissio n feature 
is a continuum subtraction (hereafter COS) method (e.g., iBouwman et al.ll200ll ). In this 



approach the continuum is fitted by a polynomial (usually first or second order) outside 
the feature and then subtracted from the spectrum. In the 5-36 yum domain {Spitzer IRS) 
there are three possible wavelength ranges in which the continuum can be fitted; shortward 
of 7.5 /im, between 12 and 15 /im and longward of 30 /im. The wavelength domains for 
the continuum fitting are not strictly determined. In this study we chose the following 
domains for the continuum fitting 6.8-7.5 /im, 12.5-1 3.5 /xm and 30-36 / xm. A reasonable 
way to subtract the continuum has been developed by Ivan Boekel et al.l (l2005bl ). Here the 
subtracted and normalized flux is given by 



^norm 



1 + 



F„ - F: 



cont 



^ Fcont \ 



Here, F^ is the flux density of the measured spectrum, F^°^^ is the flux density of the fitted 
continuum and < F^^"^^ > is the average value of the continuum flux density over all the 
frequencies in the fitted domain. For the spectral decomposition problem one can write 



N M 

EEC 

i=i j=i 



norm 



(2) 



where N is the number of dust species and M is the number of different grain sizes, nfj^"^ is 



the normalized mass absorption coefficient of the dust species i and grain size j. The mass 
fraction of a specific dust species can be calculated from the Cij values as 



a 



m. 



i=i 2^j=i 



M 



(3) 



The mass absorption coefficients should be normalized exactly the same way as it has been 
performed for the observed spectrum. In this fitting approach Kij is given, while Cij has to 
be fitted. 
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A further step was taken by adopting a single black body function as the co ntinuum 



In this single-temperature (hereafter ST) approach (e.g.. Ivan Boekel et al.ll2005bl ) the con- 
tinuum and spectral features are fitted simultaneously and the observed flux density at a 
specific frequency is given by 

N M 

= CoB^iT) + C^,JB,{T)K,,,. (4) 

i=i j=i 

Here Kij is the mass absorption coefficient of the dust species j and the grain size i, T is 
the temperature in the continuum and in the feature. B^{T) is the Planck function at a 
temperature T, the quantities N and M denote the total number of different dust species 
and grain sizes, respectively. The mass fraction of a specific dust species is given by Eq.lHl 
In this spectral decomposition problem is given, while T, Cq and Cij have to be fitted. 
An obvious limitation of this method is its simplistic treatment of the temperature structure 
of the radiating dust. 

The subsequent evolution in data quality allowed the adoption of more realistic models 
to fit the 10 /im silicate feature. By fitting the continuum and the feature temperature in 
Eq.m separately one can obtain a more proper analyzing method. In this two-temperature 



(hereafter TT) approach (IBouwman et al.ll2008l ) the flux density is given by 



N M 

i=i j=i 

Here Tc is the temperature in the continuum and Ty is the temperature in the feature. In 
all the COS, ST and TT met hods, we used a generaliz ed version of the non-negative least- 



square optimization routine of iLawson &: Hanson! (jl974j ). For the ST and TT methods a grid 



of temperatures was created first. Then the optimization routine evaluated the best fitting 
chemical composition for each temperature. During the tests the temperature step size was 
10 K in the feature and 15 K in the continuum. The minimum and maximum temperatures 
in the continuum was 20 K and 4000 K, respectively. For the feature we used 20 K and 
1500 K as the lowest and highest temperature, respectively. We assumed that dust grains 
evaporate at higher temperatures, while the lowest temperature was set to the temperature 
of the interstellar dust. In order to characterize the quality of the fit, we calculated the 
value for each temperature, with 

1 f ZTTnodcl TT'obscrvcdN ^ 

v^= ^ ySEn zl^ L (6) 

Here N^, is the number of frequencies or wavelengths and / is the number of fitted parameters, 
^observed jg ^^le flux dcuslty Calculated by the 2D RT code, F^°'^''^ is the fitted flux density 
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by the spectral decomposition method, and _pj™'' is the uncertainty in i^^observed ^^^^ ^^^q 
uncertainty of the 2D RT calculation. We simulated observations by calculating synthetic 
spectra and we used simulated observational flux uncertainties by F^"""^ = 0.001 ■ iT'observed^ 
Then the best fit chemical composition was the model, which had the lowest value. Except 
for Sec. 13. 51 we did not add noise to the spectra, therefore the fitted dust parameters do not 
depend on the value 



terror 



The basic limitation of all the above described methods is the simple handling of the 
continuum and temperature of the radiating dust. Studies of protoplanetary disks have 
shown, that the silicate emission features arise from the u ppermost optically thin disk layer 



( IChiang fc Goldreichll 19971 : iMenshchikov &: Henning||l997l ) and the silicate emission zone can 



extend from less than 1 AU to a few tens of AUs in Herbig Ae systems ( Ivan Boekel et al. 



2005al ). Since the temperature in such wide range of radii should not be constant, it seems 
to be reasonable to apply a distribution of temperatures to model the 10 fim silicate feature 
and the underlying continuum. 

Here, we propose such a new method, which is more sophisticated compared to the 
above described models, but simple enough to be applied for large data sets. In this method 
we use a simple two-layer temperature distribution (hereafter TLTD method) in order to 
describe the temperature structure of the disk. In this TLTD approach the flux density at 
a specific frequency is given by 



N M Rout 

Fu = -Pl/,cont + ^ ] ^ ^ / n -^^^(-^al 

i=l 7=1 ^ 



N M ^Rout 

//n I' 

'!,tm{r))dr, (7) 



where 

32 />Rrim 



,cont = Co-^B,{T,)+cJ B,{T,Ur))dr + 

" JKin ^ 

/.Rout 

C2 / —B,{T^,^{r))dr. (8) 

J Rrim ^ 



Here i?j„ is the inner radius of the disk, Rrim is the outer radius of the puffed-up inner rim. 
Rout is the outer radius of the disk, while r is the distance to the central star and i?^ is the 
radius of the star. is the effective temperature of the central star, while Trim{r), Tatmir) 
and Tmid{f) are the temperature of the inner rim, disk atmosphere and the disk midplane, 
respectively. The mass fraction of a specific dust species is given by Eq.[3l The temperatures 
in the disk are given by 

-^rim('") -^rim,max I 1 (9) 

\ -Ttin / 



- 10 - 



-^atm('") -^atm,max ( j-, ) (■^^) 



rmid(?^) = Tmid.max ^ (H) 

V -KRim / 

In this fitting approach R^, are given, while Tnm.max, T'mid.max, Tltm.max, grim, gmid, gatm,Co, 
Ci, C2 and Cij have to be fitted. Apart from these parameters we also need to account for 
three radii (-Rin, -Rrim, and -Rout)- 

Although there are three radii there are only two free parameters, namely their ratios 
i?out/-Rrim and i?rim/-Rin- The rcasou for that is the following. The replacement of boundaries 
in Eq.[7]and Eq.[8]from R^^ to a ■ -Rin and from -Rout to a ■ -Rout will not change the resulting 
dust composition. This behaviour is due to the Ci and Cij factors, which become ■ Ci 
and • Cij in the above described case. Since the mass fractions of different dust species 
is defined by Eq.|3], rrii j does not depend on a, although the total dust mass does. This 
means that, in a certain sense, the actual values of -Rin and -Rout are irrelevant, and the 
important parameter is rather the ratio between them. Since we do not apply radiative 
transfer calculations in the TLTD method, the radius and the temperature are not coupled 
to each other in a self consistent way and the role of the boundaries in Eq.[7] and Eq.[8] is 
only the weighting of different temperatures. 

Thus, one can constrain -Rout/-Rrim and Rnm/ R\n using the reasonable assumption that 
the outer radius of the disk is larger than the outer radius of the silicate emitting region. In 
this case, for each value of gatm one can calculate the -RoutZ-Rrim where the contribution of 
the outermost annulus to the total flux at 10 /xm equals to 0.1 % assuming a constant value 
for -Rrim (the same holds for gnni and -Rrim/-Rin)- One can then neglect the flux contribution 
of the annuli behind this radius. 

This can be the most easily done by calculating the integrals in Eq.[7H8] over tempera- 
tures instead of radii. Thus, we can rewrite Eq.[7]in the form of 



Latm.mm 2-/]- 2-qatm 



X — ^ X — ^ / /I ^ — qatm 

Fu = F.,cont + ^^'J'^^'J- L -^B,{T)T^^dT, (12) 

i=i j=i 



at m, max 



where 

A:. = ^ fi „ (13) 

^ atm,max ^'m 

The value of -Rin is not important here, its role is only a dimensional matching, therefore, we 
can choose it arbitrarily, for instance to 1 AU. Eq.[8]can also be written in the same manner. 
The upper limits of the temperatures (Trim.max, ^atm.max and Tmid.max) are fitted. Using the 
above described assumption, one can calculate the minimum temperature for each component 
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(^rim.min, T'atm.min and Tmid,min), wheie the Contribution of the corresponding annulus to the 
total flux at 10 /xm equals to 0.1%. The mass fraction of a given dust species is given by 
niij = Dij/ ^i^i ^iJ- order to compare the temperature distribution estimated by 
the TLTD method to the calculated distribution with the 2D RT code, we used Eq.[7] and 
Eq.[8] during the tests and the same i?in and i?out values were used in the fltting routine as 
in the 2D RT code. 

The fitting of the temperatures at the boundaries (Tj) and the power law indices (gj) 
is a highly non-linear problem which is not straightforward to solve. Beside the non-linear 
nature of the problem, the surface can have numerous local mi nima, which make s the fit 
even more difficult. Most of the fitting routines like, e.g., amoeba (jPress et al.ll200ll ) are not 
able to find the global minimum or/and are very ser isitive to the initial values. Therefore, 
we used the genetic optimization algorithm pikaia ( ICharbonneaul Il995l ) in order to fit the 
temperatures and the power indices. It has been shown that this algorithm can efficiently 
find the global minimum ev en in a high dimensional parameter space with multiple local 
minima (ICharbonneaul Il995l ) . This algorithm performs the maximization of a user defined 
function, thus we used the g{T,q) = l/x^ function for that purpose. Once the temperature 
distribution is fixed, the evaluation of the Cij values is a linear p roblem. Thus, we used fo r 
that purpose the generalized linear least-square fitting routine of iLawson &: HansonI (jl974j ). 



3. Results 



In the following, we investigate the quality of the decomposition methods. We fitted 
our 2D RT model spectra between 7 and 14 /xm since this is a widely used wavelength range 
for spectral decomposition in the case of Spitzer observations. In order to measure the 
quality of the fits we calculated the mass-averaged grain size of the amorphous material and 
the crystallinity fraction from the chemical composition estimated by the fitting methods. 
We use these parameters, since grain growth and crystallisation are among the key processes 
occurring during the evolution of protoplanetary dust. The evolutionary status of dust grains 
is, therefore, usually characterized by these two parameters which are commonly derived from 
the observations. For a more straightforward analysis of the results, we defined a normalized 
crystallinity as the fraction of silicate mass in crystals. The total crystallinity is the fraction 
of total dust mass in crystalline silicates. The advantage of the normalized crystallinity is 
that, for a given, fixed dust composition, its value does not change if we exclude carbon from 
the fitting procedure. 
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3.1. Disk structure and stellar type 

In the following, we investigate the effect of the disk structure and the stellar type on 
the fitted dust composition. There are three questions we address in this section, (i) How 
large is the average deviation of the dust parameters fitted by a specific method from the 
input parameters? (ii) How large is the scatter around the mean dust parameters? (iii) Is 
there any correlation, within the scatter, between the dust composition and the flaring index, 
the inclination angle or the stellar type? In order to quantify the answers for these questions 
we compute the mean value, the standard deviation around the mean and the Pearson's 
correlation coefficient between the chosen parameters and flaring index, inclination angle 
and the stellar temperature. The Pearson's correlation coefficient (r), which is frequently 
used to measrue the linear dependence of two variables {x,y), is calculated as 

^ -^iUi ■^i "^^i Vi (14) 



xf - {Y, Xi)^\/n yf - (E Vi 



where n is the number of data points. In order to investigate the significance of the Pearson's 
correlation coefficient, one can calculate the probability that the observed relation can be 
produced by random distribution with the same sample size (n) 

p{r, n) = J:lJ / (1 - uy^-^y^du (15) 

[ — ) J\r\ 



see, e.g.. lTaylorlll997l ) 



The normalized crystallinity values obtained by the COS method show a large scatter, 
overestimating and (in a few cases, e.g., for Brown Dwarf in CS series) underestimating 
the real value (see Fig.[2]Le/t). On average, the normalized crystallinity is overestimated 
by a factor of 2.8±1.9 and it correlates with the stellar temperature (r=0.53, p=10~^) and 
with the flaring index (r=0.47, p=5 ■ 10~^). The ST method overestimates the normalized 
crystallinity (by a factor of 4±1.1), which does not depend on the stellar temperature, but 
it correlates directly with the flaring index (r=0.63, p=3 ■ 10~^). There is also a hint of a 
weak inverse correlation between the crystallinity and the inclination angle (r=-0.3, p=0.01). 
The TT method overestimates the normalized crystallinity by a factor of 1.8±0.25, which 
depends weakly on the stellar temperature (r=-0.34, p=0.003) but is independent of the 
flaring index or of the inclination angle. 

The total and normalized crystallinity fraction estimated by the TLTD method is the 
closest to the input value in general. This method overestimates the crystallinity by a factor 
of 1.2±0.13. The crystallinity does not depend on the stellar type but it depends weakly 
on the flaring index (r=-0.47, p=2 ■ 10^^). The effect of the inclination on the fitted dust 
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composition is the weakest for the TLTD method. Within the studied range of inchnation 
angles (0°-70°) the scattering in the crystalhnity due to the inchnation angle is less than 1 % 
for the TLTD method. It can also be seen in Fig.^Left that the inclination has the smallest 
effect on the results of the TLTD method, indicated by the error bars, compared to other 
methods. 

Although, the input dust composition to the 2D RT disk models did not contain large 
(6 fim) enstatite grains at all, the bulk of the crystalline silicates, predicted by the TT 
method (~60 % in terms of the total crystal mass), is in the form of large enstatite grains. 
The TLTD method also suffers from the overestimation of the large enstatite fraction, while 
the ST method has this problem just in a few cases. The overestimation of large enstatite 
content is probably not related to the type of the decomposition method, since all the 
methods have this problem, but it is probably related to the similarity between the optical 
data of large enstatite and other dust species. The mass absorption coefficients of the large 
enstatite grains can be reproduced, in the fitted wavelength range, by a linear combination 
of small olivine, small carbon and medium sized enstatite grains. If this problem is really 
caused by the degeneracy of the optical data, we should get better results if we use a broader 
wavelength interval for the fits. The probability, that two dust species have similar optical 
data, is lower for a broader wavelength range. We will discuss the effect of the fitting range 
in more detail in Sec. 13. 31 

The mass-averaged grain size of the amorphous silicates has also been calculated from 
the fit results {Fig.^Left) . Since this quantity can change by more than an order of magnitude 
(0. 1-6.0 yum) we used the standard deviation of the logarithm of grain size instead of the 
grain size itself to quantify the average and scatter of this parameter. The COS method 
overestimates this quantity in all cases by a factor of 2.6^^''j'. This fitting approach predicts 
larger grains for lower stellar temperatures (r=-0.53, p=5 ■ 10^^) and for fiatter disks (r=- 
0.46, p=5 • 10"^), but the results do not depend on the studied range of inclination angles 
{0° < i < 70°). The ST method underestimates the mass-averaged grain size by a factor of 
0A2^qI. The results of the ST method depend only on the stellar temperature, predicting 
larger grains for lower stellar temperatures (r=-0.51, p=4 ■ 10~^). The TT method also 
underestimates the average grain size for T Tauri and Brown Dwarf disk models, while 
the fitted grain size is larger than the input value for Herbig Ae spectra. On average the 
TT method underestimates the average grain size by a factor of 0.7 ^'^05, showing a larger 
scatter around the mean value and correlating weakly with the stellar temperature (r=0.51, 
p=4 ■ 10~^). The TLTD method overestimates the mass-average grain size of the amorphous 
material by a factor of 1.8 and the grain size also correlates with the fiaring index (r=- 
0.6, p=10~*). It is interesting to note, that this correlation decreases towards lower stellar 
temperatures (r=-0.7, r=-0.6 and r=-0.3 for Herbig Ae, T Tauri and Brown Dwarf models. 
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respectively) . 



3.2. Featureless grains 



The amount of "featureless" grains (like, e.g., carbon) can also be important in evalu- 
ating the fraction of crystalline grains and mass-averaged grain size. Due to their smooth, 
featureless opacity curve in the fitted region (7-14 yum), the optically thin emission of such 
grains acts like a continuum. Thus, there is a degeneracy between the optically thin emis- 
sion of featureless grains and the optically thick emission of different parts of the disk. This 
degeneracy can have a strong effect on the results of the fits. The importance of this de- 
generacy is well reflected by the fact, that the ST method estimated a very high fraction 
of amorphous carbon (up to 90 % in terms of total dust mass) in all test cases. The TT 
method also had this problem in several cases, while the TLTD method did not show such 
a behavior. We discuss this behaviour in more detail in Sec. HI 

In order to study the role of carbon grains in our fits, we tried to fit the same model 
spectra without carbon (see Fig. [3]), however carbon is present in the input dust composition 



1j. We want to note, that real observations were always fitted without carbon if the COS, 
the ST or the TT method has been used for spectral decomposition. In the case of the COS 
method the resulting crystallinity is completely underestimated compared to the input value 
(down to % for Brown Dwarf spectra, see Fig.^Left) . Moreover, the COS method fits lower 
crystallinities for lower stellar luminosities (r=0.75, p=10""^^). The resulting crystallinity 
fraction of the ST is overestimated by a factor of 4±1 compared to the 'real' value. The 
resulting crystallinity in this case depends weakly on the stellar temperature (r=-0.4, p=6 • 
10~^) and on the fiaring index (r=-0.43, p=1.5 ■ 10~^) and the fitted crystallinity is higher 
compared to the fits with carbon. The TT method underestimates the crystallinity slightly 
(by a factor of 0.85±0.2). The results are sensitive to the inclination of the disk, indicated 
by the larger error bars in Fig.^Left. The results do not correlate with the fiaring index, or 
stellar temperature but depend weakly on the inclination angle (r=0.43, p=1.5 ■ 10~^). The 
average crystallinity evaluated by the TLTD method is overestimated by a factor of 1.5±0.4 
on average. The normalized crystallinity fitted by the TLTD method is less sensitive to the 
inclusion or exclusion of featureless grains compared to other methods. 

The mass-averaged grain size of the amorphous silicates is presented in Fig.^Right for 
fits without carbon among the fitted dust species. The exclusion of the featureless dust 
species results in higher average grain sizes for the silicate species for the COS and TT 



In reality carbon grains can very well be an important grain component of protoplanetary disks. 
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methods by a factor of 2 and 1.4, respectively, compared to the fits with carbon. There is a 
strong inverse correlation between the average grain size, predicted by the TT method, and 
the stellar temperature (r=-0.42, p=2 ■ 10~^) and the flaring index (r=-0.66, p=2 ■ 10"^'^). 
In the case of the ST method the resulting mass-averaged grain size of the amorphous 
silicates decreases slightly (by 20 %) compared to the fits with carbon. The average grain 
size estimated by the TLTD method decreases on average compared to fits with carbon by 
20 %. Furthermore, the the predicted grain size correlates with the flaring index (r=-0.49, 
p=8- 10-6). 



3.3. Wavelength range 

The wavelength range, for ground-based observations, is usually limited to 8-13 fim 
due to the strong absorption of the Earth's atmosphere both shortward and longward of 
this range. For Spitzer data, the broadest possible wavelength range is determined by the 
instrumental capabilities of the IRS instrument (5.2-38 /im). It may be important how 
the wavelength range for the fit is chosen. A narrower wavelength range corresponds to a 
narrower range of temperatures one has to take into account in the modelling of the featur^. 
On the other hand, a narrower wavelength range can result in a higher degeneracy among 
the optical data of different dust species. This effect results in a higher uncertainty in the 
derived dust parameters. In order to investigate how the results depend on the choice of 
the selected wavelength range, we fitted all the T Tauri spectra with different wavelength 
domains. Three wavelength regions were used for the ST and TT methods (8-13 /im, 7- 
14 /im, 7-17 /im), while these ranges were supplemented by a fourth very broad range for the 
TLTD method (5-35 /xm). 

The results showed that, for the COS, ST and TT method, the 7-14 /im range gives the 
best results in general. However, the estimated large enstatite particle fraction is much closer 
to the "real" values, using a broader wavelength range (see Fig.^iight) . The estimated mass- 
averaged grain size of the amorphous silicates is underestimated with the other wavelength 
intervals (roughly by a factor of two) compared to the results of the fit using the 7-14 fim 
domain (see Fig.^iight). In the case of the TLTD method, the broader wavelength range 
results in better estimates of all parameters. It is interesting to note that we get the best 
results using the 7-17 /im range and the results become worse if we use the broadest range 
between 5 and 35 yum. Using the 7-17 /zm wavelength interval the mass-average grain size 



^For a narrower temperature range, the assumption of a fixed dust composition is certainly better justified 
compared to a wide range in temperatures which corresponds to a larger disk region. 
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and total crystallinity, fitted by the TLTD method, are overestimated by by a factor of 1.45 
and 1.05, respectively. In this case, the only correlation we found is the dependence of the 
mass-averaged grain size on the flaring index of the disk (r=-0.7, p=8 ■ 10~^^) Using the 
broadest wavelength range, the average grain size is underestimated by 70 %, while the total 
crystallinity is overestimated by a factor of 1.75. One reason for that can be that even the 
TLTD method is too simple to describe the temperature structure of a real disk correctly. 
Another reason for the worse results can be found in the optical data of the different silicate 
species. While the crystalline silicates have sharp features in the 20-30 fim region, the optical 
data of large amorphous silicate grains are very smooth, which makes it difficult to recognize 
this component. 

The normalized crystallinity as a function of mass-averaged grain size is presented in 
Fig. for the best configuration of each method. The results of the TT method seem to be 
relatively close to that of the TLTD method. However, the fitted dust composition obtained 
by the TT method suffers from the overestimation of the abundance of large enstatite grains. 
The predicted amount of large enstatite grains was closer to the input value using a broader 
wavelength range (7-17 fim) regardless of the spectral decomposition method. This fact 
proved our former assumption that the overestimation of the large enstatite grains is the 
result of the degeneracy between the optical data of the large enstatite grains and other dust 
species in the narrower wavelength interval. 



3.4. Large grains 

It is well known, that the strength of the spectral features in the 10 /im region decreases 
with increasing grain size (see, Fig.[T]). Therefore the emission of dust grains grains larger 
than 2 fim are usually neglected during the spectral decomposition. Although the strength 
of the 10 /im feature of the amorphous silicates is much weaker for a grain size of 6 /xm than 
for 0.1 /im the feature strength is not completely zero. In this study, we included also 6.0 fim 
sized grains in the input dust composition to the 2D RT code. Dust grains with such size 
are probably present in the disk atmosphere if grain growth occurs in the disk, although 
their abundances depend on the turbulence and sedimentation processes. In order to test 
the importance of the largest grain population (in case of ground-based observations), all 
T Tauri spectra were fitted using a wavelength range of 8-13 fim and excluding all 6 fim 
sized grains from the fit, although large grains were present in the 2D RT disk models. In 
this case we were interested how well the properties of the smallest grain populations can 
be recovered by the methods. Therefore, during the comparison of the fitted dust mixtures 
with the input composition we used only the 0.1 fim and 1.5 fim sized grains, which is the 
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usual approach in the hterature. 

The results are summarized in Fig.O It can be easily seen that the dust parameters 
fitted by all methods deviates from the input dust parameters. If we included large grains in 
the fit, their estimated fraction by the ST method was less than a few percent, therefore their 
effect on the fit is obviously weak. The resulting dust parameters of the TT method, get 
further away from the right values if we excluded the large grains from the fit. The smaller 
mass-averaged grain size is caused by the excluded large amorphous grains, the fraction of 
which is significant if we included them in the fit. The smaller estimated crystallinity is a 
result of the exclusion of large enstatite grains. The results of the TLTD method were almost 
identical to the input dust mixture giving far the best results among the tested methods. The 
differences between the resulting dust composition including/excluding large grains in/from 
the fit can be related to the disk structure and will be discussed in detail in Sec. 14. 31 



3.5. The influence of simulated noise 

Observations always suffer from uncertainty in the observed quantities and these errors 
naturally affect the parameters derived from the observational data. Since one cannot avoid 
the effect of noise, it is important to know how sensitive the applied method is to the 
noise level and how the uncertainty of the derived parameters can be estimated from the 
observational errors. In the case of the spectral decomposition of the silicate feature around 
10 /im, it is very difficult to describe the error propagation by analytical formulas. A simpler 
way to estimate the u ncertainty in the fitted c hemical composition is a Monte Carlo type of 



error estimation (e.g., Ivan Boekel et al.ll2005bl ). In this kind of error estimation, a normally 



distributed noise is added to the spectrum, scaling the width of the distribution to the 
simulated observational uncertainty in the flux value. Then the resulting spectrum should 
be fitted. This procedure should be repeated many times. Then the standard deviation of 
the resulting mass fractions will be the uncertainty of the derived dust compositions. 

We used the spectrum of a T Tauri star with a moderately flared disk and an inclination 
of 45° for testing the effect of different noise levels on the results of the TLTD method. We 
studied the effect of noise on ground- and space-based observations, using a wavelength range 
of 8-13 /im and 7-17/im, respectively. Three different noise levels have been used in these 
tests with an i?error/ ^observed q q ^^^^ q qq^^ Corresponding to a signal-to-noise (S/N) 

ratio of 10, 100 and 1000, respectively. We want to note, that the noise level was assumed to 
be independent of wavelength, which can be good approximation for bright sources, while this 
assumption is certainly not correct for faint sources with strong silicate feature. After 100 
spectra have been generated and fitted for each noise level, we studied the average chemical 
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composition over the 100 spectra, the mass-averaged grain size and normahzed crystalhnity 
and the scattering around the average values. The latter two parameters have the advantage 
that, even for low S/N, they can be more certainly evaluated than the mass fraction of a 
single dust component. 

The normalized crystalhnity as a function of mass-averaged grain size, derived by the 
TLTD method, is presented in Fig. [7] for simulated ground-based observations (8-13 /im). It 
is easy to see , that the uncertainty in the value of crystalhnity and mass-averaged grain size 
(obviously) increases with increasing noise level. Moreover, for lower signal-to-noise ratio 
the average value of the crystalhnity and the average grain size increases. For a S/N of 100 
the average value of the crystalhnity and the mass-averaged grain size are overestimated 
by a factor of 1.5 and 3, respectively, compared to the input values. The scattering of the 
individual points is not perfectly symmetric to the mean value. The uncertainty of the value 
of the crystalhnity is higher above the mean value than below that. The uncertainty on the 
value of crystalhnity for a S/N of 100 is +11% -7% in an absolute sense (i.e. in terms of 
the total silicate mass). In the case of the average grain size the uncertainty on the average 
value is +1 /im and -1.5 /im. 

For space-based observations (see Fig.[8|, 7-17 /xm range), the scattering of the normal- 
ized crystalhnity and mass-averaged grain size decreases compared to that of the ground- 
based observations. Furthermore the average values of these quantities are closer to the 
input value of the 2D RT codes using the broader wavelength range. For a S/N of 100 the 
average value of the crystalhnity and the mass-averaged grain size are higher than the input 
value by a factor of 1.26 and 1.44, respectively. The uncertainties on the average value of 
the crystalhnity is 5 % in an absolute sense (i.e. in terms of the total silicate mass). In the 
case of the average grain size the uncertainty on the average value is 0.8 fim. 

Although the signal-to-noise ratio of the data is of fundamental importance, the un- 
certainty of the fitted dust parameters, caused by the noise, depends on the actual dust 
composition and the optical constants. Since the crystalline silicates have usually strong, 
narrow and sharp features their abundances can be constrained better, even for low signal- 
to-noise ratio, than the abundances of the amorphous silicates with broad, flat features. 
Therefore, it is obvious, that the higher the crystalhnity the easier it is to determine its 
value. 
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4. Discussion 

Our results clearly indicate that the evaluation of the chemical composition of dust in 
protoplanetary disks from mid-infrared spectra is not straightforward even if do not have 
to consider the uncertainties in the optical constants of the dust grains. The resulting dust 
parameters strongly depend on the assumed continuum below the feature and the assumed 
temperature (s) of dust grains emitting at the considered wavelengths. In the following, we 
will discuss the reasons why the results of different methods deviate from the input dust 
composition, and which effects should be taken into account during the analysis of mid- 
infrared dust features. 



4.1. Temperatures in the feature 

Both the ST and TT methods assume that the source function underlying the emission 
in the 10 feature can be well app roximated by a blackbody curve with a single effective 



temperature. Ivan Boekel et al.l (j2005al ) have shown that the region where the 10 /xm emission 
comes from can extend from less than lAU to several tens of AU in Herbig Ae systems. 
They also have shown that the mid-infrared radiation, emitted by the different annuli (with 
different temperatures), contributes to the total observed flux in this wavelength regime 
significantly. We also investigated the temperature structure of the 2D RT disk models 
in the region where the bulk of the 10 /xm emission originates. Images of the disk models 
were calculated at 10 fim and then the cumulative flux as a function of radius was used to 
determine the size of the region from which 70% of the total flux originates (Fig.lHh)- In 
order to characterize the temperature in the disk atmosphere we evaluated the temperature 
along the li ne where the radial optical depth ro.ss^m reaches unity in the disk, using the 



definition of IChiang fc GoldreichI (119971 ). From Fig.[n]6 it can be seen, that the temperature 
in the disk atmosphere changes significantly in the region where the silicate emission comes 
from. On the other hand the changes of the surface density of the disk atmosphere as 
a function of radius are just marginal and they depend on the flaring index. Since the 
silicate emission feature is produced by a convolution of optical depth (or surface density) 
and the Planck function at different temperatures, it now becomes obvious why the simple 
"isothermal" models give poor results. A sum of different Planck curves with these rather 
different temperatures cannot be reproduced by a single Planck curve. 

Since we used the same Rm and Rout both in the 2D RT code and in the TLTD fitting 
method, we can compare the fitted temperature distribution with the temperature structure 
of the 2D RT models. Fig.[ni> shows that the TLTD method can efficiently evaluate the tem- 
perature distribution in the disk atmosphere. The ratio of the slope of the fitted atmosphere 
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temperature by the TLTD method and that of the 2D RT code was 0.88±0.07 averaging over 
all the fitted spectra. The same ratio for the boundary temperature (Tatm.max) was 1.3±0.6. 



4.2. Temperatures in the optically thick continuum 

Another assumption of the ST and TT methods is that the continuum below the emission 
feature can be described by a Planck function w ith a single temperature. Using the passive 



irradiated disk model of iDuUemond et al.l (120011 ) . the continuum emission in that region has 



four components: emission of the central star, the puffed-up inner rim, the disk midplane 
and the emission of the featureless (e.g., carbon) grains in the disk atmosphere. In general, 
shortward of ~8 /im, the continuum is dominated by the emission of the hot inner rim of 
the disk with decreasing flux towards longer wavelengths. At wavelengths longer than 12- 
14 //m most of the continuum emission comes from the disk midplane with increasing flux 
towards longer wavelengths. Thus, the slope of the continuum should change in the 8-12 /xm 
region in order to match the decreasing flux of the inner rim and the rising flux of the 
midplane. The turn-over point, where the slope changes is determined by the balance of the 
rim and midplane components (which is set by the disk geometry and the inclination) and 
by the featureless grain (e.g., carbon) content of the dust. One of the reasons for the poor 
performance of the COS, ST and TT methods lies in their incapability to fit a continuum 
with such a shape. 

The ST method simply cannot reproduce the real continuum, since the temperature in 
the continuum and in the feature is the same in this approach (see Fig.lTOb). The only way 
to fit the real continuum is the inclusion of large amorphous carbon, the mass absorption 
coefficient of which increases with wavelength in this region. Therefore, this method fitted an 
unrealistic amount of carbon (~70-90 % of the total mass) in order to fit the real continuum. 

In the case of the TT approach two distinct solutions exist for the fit, which can be easily 
seen in Fig.[TT]as the space has two minima. This kind of space is characteristic for the 
TT metho d during the fitt i ng of the spectrum of a Class II source, using the classification 



scheme of iLada fc WilkingI (119841 ) and related to the fit of the continuum. One of the two 
minima has a high continuum temperature (~ 600 K) and low feature temperature (~ 200 K, 
see Fig.fTOb). In this case, the resulting dust composition overestimates the mass fraction of 
large (a = 6 ^m) carbon grains. In this solution the real continuum is fitted with the high 
temperature continuum component shortward of ~10 fim, while at longer wavelengths the 
real continuum is fitted by the optically thin emission of carbon grains. In the case of the 
other minimum, the feature has a higher temperature than the continuum and the resulting 
amount of carbon is very well estimated. On the other hand the silicate grains, producing the 
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10 jjm feature, are significantly larger compared to those in the other solution. The reason 
for the larger average grain size is that the mass absorption coefficient of the large (6 fim) 
amorphous grains is higher than that of the smaller ones shortward of ~8.5 fim. Therefore the 
real continuum emission at longer wavelength is fitted with the single temperature continuum 
of the TT method, but shortward of 8.5 /im the real continuum is fitted by including large 
amorphous grains. 

The subtraction of the continuum can be a better solution if, and only if, the real 
continuum is known, otherwise the subtraction of the continuum introduces an additional 
uncertainty in the results. The highest uncertainty in the COS method is that it tries to 
extrapolate the continuum by fitting a few data points. Therefore, the fitted continuum 
is very sensitive to the wavelength ranges in which it has been fitted. Although, the COS 
method estimated the mass-averaged grain size in a few cases equally well as the TLTD 
method, the total crystallinity was deviating most strongly from the real value. 

We compared the temperature profile of the disk midplane calculated by the 2D RT 
code to the fitted temperature distribution by the TLTD method (Fig.lHlc). In this case 
the temperature distribution in the 2D RT model was calculated along the line where the 
vertical optical depth in the disk at 10 fim equals unity. One can see that the TLTD method 
overestimated the temperature compared to that of the 2D RT code. Furthermore the 
slope of the temperature distribution is very low. The ratio between the fitted power-law 
indices and the calculated one by the 2D RT code was 0.33±0.33, and that of the boundary 
temperature was 1.4±1.0. The reason for the difference between the fitted and the 'real' 
midplane temperature in general is caused by the fat that even the TLTD method is too 
simple to represent the temperature structure of a protoplanetary disk. In contrast to the 7- 
8 fim region, longward of 12 fim the continuum cannot be unambiguously determined, since 
the dominant contributor to the total emission is still the disk atmosphere. Moreover, there is 
a degeneracy between the optically thin emission of the featureless grains and the optically 
thick disk emission. In a 2D RT code the layer where the 'midplane' emission originates 
depends strongly on the disk geometry, inclination angle, and of course on the wavelength. 
During the comparison of the results of the TLTD method and the 2D RT code we made the 
simplifying assumption that the layer, where the 'midplane' emission comes from, is where 
the vertical optical depth at 10 fim equals unity, which is only true for face-on orientation. 

4.3. Explanation of the correlations 

The shape and strength of the silicate feature around 10 /im is affected by the disk 
structure, the inclination and the stellar type in a complex way. The efficiency of a spectral 
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decomposition method depends on how it can handle these effects. In the former two sections 
we discussed why the single and two temperature approximations or a continuum subtrac- 
tion are insufficient to model the temperature structure of the silicate emission zone in a 
protoplanetary disk. This weakness of the fitting routines can lead to a large scatter in the 
derived dust parameters. Here, we discuss how and why the derived dust parameters can to 
depend on the geometry of the source or the stellar type. If such an "artificial" dependence 
is not properly understood, the comparison of grain parameters derived for different source 
types can be misleading. 

The effect of the stellar temperature on the mid-infrared silicate features can be much 
stronger than it is usually thought. Under the assumption of an identical star-to-disk lu- 
minosity ratio and disk structure, the contribution of a Brown Dwarf (Tcff=2500 K) to the 
10 /xm emission is roughly 44 times higher than that of a Herbig Ae star (Teff=9500 K). In 
Fig.[T2la it can be seen that this increasing "stellar" contribution toward lower effective tem- 
peratures results in lower peak over continuum ratios of the 10 /im feature without changing 
the parameters of the dust model. This phenomenon naturally explains the strong inverse 
correlation between stellar temperature and the temperature fitted by the ST method (r=- 
0.84) and the continuum temperature fitted by the TT method (r=-0.92). Moreover, the 
feature temperature over continuum temperature ratio predicted by the TT method directly 
correlates with the stellar temperature (r=0.93). The continuum temperature is always 
higher than the feature temperature for Brown Dwarfs and mostly for T Tauri stars. This 
kind of soluti on seems to be unphysica l (compared with the radiative transfer solution of 



the problem (jChiang &: GoldreichI 119971 : iMenshchikov fc Henning 119971 )). The TT method 



is not a solution of the radiative transfer problem and simply fits the dominant continuum 
component. As the stellar contribution becomes the dominant continuum component for 
lower stellar temperatures, the TT method fits higher continuum temperatures. This be- 
haviour of the TT method leads to a strong "artificial" inverse correlation between the stellar 
temperature and the mass-averaged grain size of the amorphous silicates if we fit the spectra 
without carbon. The TLTD method is much less affected by this effect, since the stellar 
contribution is already included in the fit. The average grain size does not depend on the 
stellar temperature if we include carbon in the fit, but it does if we exclude carbon, although 
this correlation is weaker compared to that of the other methods (r=-0.67). 

The effect of the flaring i ndex, i.e. the disk geomet ry, is more complex than that 



of the stellar temperature. As iKessler-Silacci et al.l (120071 ) have shown the radius of the 



zone where the 10 /zm emission comes from can decrease for flattened disks compared to 
flared disks. Thus the decrease of the flaring index can have two important consequences 
concerning the 10 /im silicate emission band. One effect is that the contribution of the 
midplane component to the total flux below the silicate feature decreases toward flattened 
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disks ( iDullemond fc Dominikll2004l ). In other words, the dominant continuum temperature 
shifts toward higher values toward lower flaring indices. Moreover, the amount of opti- 
cally thin matter decreases more rapidly with radius in flattened disks than in flared disks 



( iKessler-Silacci et al.l 120071 ) . This behaviour leads to the decrease of the height of silicate 
emission feature above the continuum toward flattened disks (see Fig.[T2l6). The 10 /xm fea- 
ture is, however, less affected by the disk geometry for Brown Dwarfs compared to Herbig 
stars, since the relative contribution of the "stellar" component to the continuum emission 
at 10 fim becomes stronger for lower stellar temperatures. This behaviour explains why the 
temperature fitted by the ST method inversely correlates with the flaring index for Herbig 
Ae spectra (r=-0.85) but not for Brown Dwarf spectra (r=-0.33). The continuum tempera- 
ture as well as the continuum over feature temperature ratio predicted by the TT method 
behaves exactly the same. This is likely the reason why the mass-averaged grain size in- 
versely correlates with the stellar temperature for all the methods (however the strength of 
the correlation depends on the method), but this correlation is the strongest for Herbig Ae 
spectra and almost disappears for Brown Dwarf spectra. The TT method is an exception 
as the predicted grain size correlates inversely with the flaring index also for Brown Dwarf 
spectra. 

The results of none of the methods show a close correlation with the inclination angle. 
The scatter in the derived dust parameters is the smallest for the TLTD method (more than 
a factor of two compared to other methods). It can be seen in Fig.[T2b) that the changes 
in the inclination angle modifies the relative contribution of different (mainly continuum) 
components (i.e., star, inner rim, and disk midplane) to the total flux. The high temperature 
continuum components (the star and the inner rim) becomes stronger for higher inclination 
angles (close to edge-on) while at the same time the midplane component becomes weaker. 
The TLTD method can better handle this behaviour due to the Cq, Ci, C2 factors. We want 
to note, that the TLTD method, in this form, cannot handle edge-on disks where the silicate 
emission band is affected by the extinction of the outer parts of the disk. 

In general, all these findings should serve as a warning, that correlations between grain 
parameters and disk/stellar properties can be artificially produced by limitations of the 
spectral decomposition methods. 



4.4. Degeneracies in the optical data 

Another problem in the spectral decomposition of the mid-infrared silicate feature can 
be the degeneracy among the optical data of the fitted dust species. If the mass absorption 
coefficients of two dust species are very similar to each other it will be hard to identify 
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them separately in the observed spectrum regardless of the applied spectral decomposition 
methoco- 



Similar to Ivan Boekel et al.l (l2005bl ). we also investigated this effect by fitting the mass 
absorption coefficients of each dust species with that of the other species. In order to quantify 
the quality of the fit we calculated the average deviation of the best fit model from the mass 
absorption coefficient of the analyzed dust species (< >), where 



model 

V 



1/2 



-^=^^^T-^T^ • (16) 



Here, is the mass absorption coefficient of the dust species under consideration, and 
is the mass absorption coefficient of the best fit model. The value for each fit is presented 
in Fig.[ 



It can be seen, that for a wavelength range of 8-13 /im the optical data of pyroxene- 
type amorphous silicates and large (6 /xm) enstatite can be well reproduced by the other 
dust species. For this wavelength range there are several dust species for which the < cr > 
value is less than 0.1, implying that they can be on average reproduced at a 10 % level by a 
linear combination of other dust species. If one has a signal-to-noise ratio of 10 these dust 
species cannot be identified in the spectrum due to the similarities of the mass absorption 
coefficients. It can be seen from Fig. [12] 6 and Fig.[T3]c that < a > increases with the width 
of the wavelength range used for the fit. On the other hand, the comparison of Fig. [13] c and 
Fig.[T3]c? shows that, for a wavelength range of 5-35 /im, the optical data of a specific dust 
species can be equally well reproduced by the other species, as for a wavelength range of 7- 
17/im. If the mass absorption coefficients of two dust species are very similar to each other, 
then they cannot be identified in the spectrum separately independently from the applied 
spectral decomposition method. We recommend to perform this kind of analysis for any new 
dust model, in order to recognize and understand the possible degeneracies. 

During this analysis we did not take into account the continuum since the shape of 
the fitted continuum depends on the applied spectral decomposition method and we were 
interested in the degeneracies among the mass absorption coefficients only. The applied 
continuum has probably the strongest effect on the featureless dust component (in our case 
on the carbon). Although the real continuum (the shape of which is the best reproduced 
by the TLTD method) is not a straight line as the mass absorption coefficient of the 6 jj^m 



•^Dust components will not behave the same in very wide wavelength regions. Radiative transfer cal- 
culations would then require to obtain dust temperatures which would allow to distinguish between the 
components. 
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sized amorphous carbon grains in the fitted domain, it can still be very difficult to identify 
such featureless dust species, especially using a narrow wavelength range (8-13 /im). The 
mass absorption coefficient curve of the 0.1 and 1.5/im sized amorphous carbon grains can 
also be hardly separated from the combination of inner rim and disk midplane emission. 
Therefore, one should be careful during the interpretation of the carbon fraction fitted by 
simple methods like those ones tested in this paper. 

In order to t est the significa nce of the detection of a given dust component an F-test 



can be used (e.g. iMin et al.ll2007l ). Synthetic observations should be created by adding sim- 
ulated noise to the observed spectrum as it has been done in Sec l3.5[ These synthetic spectra 
should be fitted with and without the dust component under consideration. The ratio of the 
two distribution (with and without the dust component) gives the F-distribution. After 
the mean value and the standard deviation (a) of the F-distribution has been calculated one 
can determine the distance of the mean (in terms of a) from unity, which gives the signifi- 
cance level for the detection of the dust component. If a component, which is significantly 
present in the spectrum, is added to the fit the should decrease, while for an insignificant 
detection the x^ does not change. One should, however, be careful with the F-test during 
the significance test, since the F-test assumes that all components are independent from 
each other, i.e. there is no degeneracy among the components. Without any knowledge 
on the degeneracies among the dust components one should perform an F-test not only for 
each dust species but also for all possible linear combination of them in order to draw the 
right conclusion from the test. Another possibility can be that before the F-test is done one 
performs the above described simple analysis of the mass absorption coefficients (fitting the 
mass absorption coefficient of a given dust species with a linear combination of others). Such 
an analysis gives an indication for which dust species an F-test should be done carefully. 

The general performance of a spectral decomposition method is determined by two 
important effects, which must be handled together and which act in opposite directions. One 
problem is the degeneracy /similarities of the optical data of the ffited dust species, which 
requires the broadest possible wavelength range. The other effect is linked to the ffited 
temperature (s). All the simplifying assumptions of the above described simple methods, as 
for instance the assumption of an average temperature, becomes more valid for a narrower 
wavelength interval (i.e., a narrower range of temperatures), compared to a broad one. The 
optimal wavelength range for a specific method is therefore determined by the balance of 
these two effects. Due to the assumption of a single temperature, all previously used methods 
can handle just a limited wavelength range (8-13 /im or 7-14 /xm), where the similarities of 
the optical data can play an important role for low S/N data. An advantage of the TLTD 
method over the other methods is, that one can extend the ffited wavelength interval by 
applying a temperature distribution. On the other hand, the TLTD method is still an 
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approximate method, which means that one cannot apply this method to an arbitrary broad 
wavelength range. An optimal wavelength range can be the interval between 7 and 17/im. 
For this interval, the differences between the optical data of the fitted dust species are 
roughly as high as for a wavelength range of 5-35 /xm, applying our dust model. Moreover, 
the 7-17 /im range is narrow enough to apply our approximate temperature determination 
and allows to make the assumption that the grain composition does not change dramatically 
with radius. 



4.5. Limitations of the TLTD method 

Although the new TLTD method returned dust parameters closest to the input dust 
mixture in most cases, his method has its limitations as well as any other method based 
on simplifying assumptions. During the derivation of Eq.[T2] we assumed that the optical 
depth in the disk atmosphere does not depend on the radial distance. In a real disk the 
optical thickness in the disk atmosphere depends on the radius. However this dependence is 
marginal and depends on the flaring index. Therefore, the assumption of constant optical 
thickness in the disk atmosphere as a function of radius seems to be a reasonable assumption. 
The above described version of the TLTD method does not contain reddening. Therefore, 
we do not suggest this form of the method for disks with extremely high inclination angles 
(edge-on orientation. Here, one should perform a correction for reddening before applying 
the TLTD method. We assumed also that the 10 /xm emitting region is smaller, than the 
outer radius of the disk. If this assumption is not valid, then one has to fit the outer radius 
as well. Despite these limitation we suggest to use the TLTD method, since the bulk of the 
young stars with an age of few Myr, where dust processing is thought to be strongest, fit 
within these criteria. 

Although we tested the spectral decomposition methods on spectra of numerous disk 
models, there are still untested cases. We did not include radial mixing or self-consistent 
sedimentation in our 2D models and we fixed the input chemical composition in order to 
investigate the effect of disk structure on the resulting dust composition. These assumptions 
can have an effect on the resulting chemical compositions of all methods. We are aware, that 
our results were derived from a limited set of input models, based on a common parametriza- 
tion. However, since we have done the analysis for a number of choices of parameters of the 
input model, we feel confident that our conclusions are robust. In this paper we presented 
an efficient method, which can be applied to analyze large data sets in order to obtain a first 
estimation of the dust composition in the 10 fim emitting region. The main improvement 
in the TLTD method is the application of a temperature distribution to describe the tem- 
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perature structure in the disk and the apphcation of multiple continuum components. Since 
these assumptions are a better description compared to the assumption of the other simpler 
methods, wc think, that the advantage of the TLTD method over the other methods may 
not be affected by taking into account radial mixing, sedimentation or different chemical 
compositions. 

5. Conclusions 

Out of the three spectral decomposition methods, tested in this study, the the two-layer 
temperature distribution (TLTD) method returned dust parameters closest to the input 
parameters for the calculation of the synthetic spectra. Compared to the results of the 
TLTD method, the fitted dust composition by the two temperature (TT) method showed 
larger deviations from the input dust composition, while the results of the single temperature 
(ST) method deviated even more. Although the mass-averaged grain size estimated by the 
continuum subtraction (COS) method was very close to the input value in a number of cases, 
the fitted crystaUinity of these models was far away from the 'real' values. 

There arc several reasons why the previously used spectral decomposition methods have 
limited capabilities in recovering the dust composition from the 10 /im silicate feature. The 
COS method tries to reconstruct a complex continuum by fitting a few data points outside 
the 10 lira silicate feature. This can modify the shape of the feature introducing additional 
uncertainties to the analysis. The ST and TT methods assume, that (i) the source function 
underlying the emission in the 10 /xm feature can be well approximated by a blackbody 
curve with a single effective temperature and (ii) the continuum below the feature can 
also be described by a Planck function with a single temperature. Additionally, there is 
a degeneracy between optically thin emission of featureless grains (e.g., carbon) and the 
optically thick emission of the circumstellar disk. Due to this degeneracy, if we included 
the amorphous carbon into the fit, then the ST method always overestimated the fraction 
of carbon in the spectra, while the TT method also had this problem in a number of cases. 
The TLTD suffered less from this problem compared to the other methods due to the more 
sophisticated estimation of the optically thick emission of the disk. 

The exclusion of carbon from the fitted dust species resulted in a better estimation of 
the total crystaUinity in the case of the TT method, although the predicted mass-averaged 
grain size of the amorphous silicates became too small. For the ST method the exclusion 
of carbon results in an overestimation of the total crystaUinity (up to ~70%) even more 
so than including carbon. Furthermore, the ST method predicted a higher crystaUinity 
for lower stellar luminosities if we excluded the featureless grains from the fit. The mass- 
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averaged grain size of the amorphous sihcates was somewhat larger for the TLTD method 
without carbon in the fitted dust species. 

All the methods overestimated the fraction of large enstatite grains using a wavelength 
range of 7-14 /xm. This behavior is related to the degeneracy of the mass absorption coef- 
ficients of large enstatite grains and small olivine, small carbon and medium-sized olivine 
grains. The overestimation of the large enstatite fraction decreased if we used a broader 
wavelength range (e.g., 7-17/xm) for the fit, where the degeneracy is less pronounced. 

Although the input dust mixture contained only 4 % of large grains (6 /im) in terms of 
the total mass mass, the exclusion of this grain population from the fit had an effect on the 
resulting dust mixture, using a wavelength range of 8-13 ;um (ground-based observations). 
The fitted dust parameters are closer to the input mixture for the COS, ST and TLTD 
methods if the 6 /xm grain population is excluded from the fit. In contrast the TT method 
returned dust parameters which differed more from the input parameters if we excluded the 
largest grain population from the fit compared to the results including the 6 /xm grains. In 
the case of the TLTD method, the exclusion of the large grain components modified the mass- 
averaged grain size of the amorphous dust species. The changes in the estimated fraction of 
the crystaUine dust species, estimated by the TLTD method, are just marginal. In the case 
of the COS, ST and TT methods the fitted crystallinity was more strongly affected by the 
exclusion of the 6 fira sized grains. The reason for that was partially that by including large 
grains in the fit, all methods overestimated the fraction of 6 /im sized enstatite grains. 

Here we should note that poorly performing analysis methods can lead to artificial cor- 
relations between stellar/disk properties and dust grain parameters. The poor performance 
of such methods in representing the disk spectral energy distribution is then compensated 
by introducing emission from non-existing grain components. 

The quality of the derived dust parameters depends very strongly on the signal-no-noise 
ratio of the observed spectra. For the applied set of optical constants and using a wavelength 
range of 8-13 /xm one can expect an uncertainty of about 11 % (in terms of the total silicate 
mass) in the value of the crystallinity for a S/N of 100. For the same signal-to-noise ration, 
but using a wavelength interval of 7-17 /^.m the uncertainty in the value of the crystallinity 
decreases to 5%. We want to note, however, that probably the higher the crystaUinity the 
easier it is to determine its value regardless of the noise level. This is caused by the sharp 
peaks of the crystalline material which are easier to identify than the broad feature of the 
amorphous silicates. 

Our tests showed that the optimal wavelength range of the fit is 7-14 iim for the ST 
and TT method assuming the mass absorption coefficients used in this paper. The width of 
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this wavelength interval is limited by the similarities in the optical data, which requires the 
broadest possible interval and the assumption of a single average temperature in the feature 
and in the continuum, which acts in the opposite direction. The optimal wavelength range 
for the TLTD method has been found to be 7-17/im. The application of this broader wave- 
length range, compared to the simpler methods, resulted in a better estimation of the dust 
composition. The 7-17/im interval is broad enough to increase the differences among the 
mass absorption coefficients. On the other hand, the temperature distribution is a consid- 
erably better estimation of the temperature structure compared to an average temperature. 
Although, the TLTD method provided a better estimation of the dust composition than the 
other methods, it has been found to be too simple to fit even broader wavelength intervals 
(e.g., the full Spitzer IRS domain, 5-35 /im). One should therefore use the TLTD method for 
a first estimation of the dust composition by fitting the 10 fim wavelength range, since a more 
realistic description of the radial temperature structure and changes in the dust composition 
is necessary to fit an even broader wavelength range. Another possibility would be to apply 
the TLTD method for the outer disk (e.g., 20-40 /zm interval) and the longer wavelength part 
of the Spitzer spectra and fit this range separately. 

On the basis of the results of our study, we propose the following recipe for the analysis 
the silicate emission features. 



1. Set of mass absorption coefficients. First, one should investigate how degenerate the 
fitting problem is, i.e., how the mass absorption coefficients of a given dust species can 
be reproduced by the linear combination of that of the other dust species. If there are 
two dust species with very similar mass absorption coefficients, then one should handle 
them together, since they cannot be distinguished in the fitting procedur^. 

2. Global shape of the SED and the disk structure. The results of our study show that 
the applied assumption for the temperature structure of the disk is of fundamental 
importance. Therefore, one should investigate the global shape of the SED of the 
analyzed source in order to make a reasonable assumption for t he disk structure. If 
the gl o bal shape of the SED can be desc ribed with a disk model of Ichiang fc Goldreich 
(119971 ). Menshchikov fc Henninel fll997l l or that of bullemond et all (1200 ih . then one 
should use the TLTD method for the analysis of the silicate feature. The TLTD method 
also has its limitation. In its presented implementation it can handle only disks with 
a continuous surface density distribution. For instance, if the disk has a large gap in 



"^In a real radiative transfer calculation this may be possible since the optical properties of different 
materials are never completely identical over a large wavelength range. 
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the region where the sihcate emission comes from, then one has to modify the current 
form of the TLTD method by excluding this temperature range from the fit. 

3. Spectral decomposition with the TLTD method. Before the TLTD method is apphed 
one has to choose the wavelength range for the fit. We propose a wavelength range of 7- 
17 /im for space-based observations (with the presented dust model) , while for ground 
based observations the wavelength domain is limited by the Earth's atmosphere. The 
continuum should not be subtracted from the spectrum, since this can introduce an 
additional uncertainty in the results, as our results have shown in the case of the 
COS method. We suggest to use the TLTD method in the form of integration over 
temperatures instead of radii in order to avoid assumptions for the outer disk radius. 

4. Interpretation. The uncertainties in the resulting dust parameters can be estimated 
by a Monte Carlo type of error estimation from the observational uncertainties. For 
spectra with low signal-to-noise ratio (few tens) the fitted mass fraction of a single 
dust component can be very uncertain. Therefore we propose to use robust statistics, 
like the mass-averaged grain size or the crystallinity, for low S/N data in order to 
characterize the dust properties. Our test shows that for a S/N of 100 one can expect 
an uncertainty in the value of the crystallinity of 11% in terms of the total silicate 
mass for ground based observations (8-13 /im). The expected uncertainty for space 
based observations (7-17 /xm) is 5% for the same S/N level. 

During the interpretation, one also should take into account the systematic errors in 
the fitting method. Using the optimal wavelength range for the fit, the TLTD method 
overestimates the mass-averaged grain size of the amorphous silicates by 50±30 %, 
while the total crystallinity is overestimated by 8±7%. The increase in the noise level, 
as well as a narrower wavelength range (8-13 /xm) results in larger mass-average grain 
sizes for all dust species and higher crystaUinity values than the real values. 

Th. Henning was supported in part by the National Science Foundation under grant 
PHY05-51164 (at KITP). We thank an anonymous referee for careful reading the paper and 
excellent suggestions for improving the text. 
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Disk models 



For the moderately flared and flat disks we use the following density distribution 




(Al) 
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where 



and 



Ri 



Hp{Ri. 

Rin. 



Ri: 



(A2) 



(A3) 



Here Hp is the pressure scale height in the disk and a is the flaring in d ex. F or the puffed-up 
inner rim we applied a similar model than that of iDuUemond et al.l (1200 ll ). However, the 
scale height of the inner rim was an adjustable parameter. It has been set in such a way, 
that the shadowed region of the disks relative to their inner radius have to be the same for 
all models in a given series. 
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Fig. 1. — Mass absorption coefficients of the dust species used in the tests. All dust species 
were used with three sizes: 0.1 /^m (solid line), 1.5 /^m (dotted line), 6.0/xm (dashed line) 
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Fig. 2. — Fraction of silicate mass in crystals {Left) and mass-averaged grain size (Right) for 
(a) Herbig Ae star, (b) T Tauri star and (c) Brown Dwarf with carbon among the fitted dust 
species. The symbols show the averaged value over the different inclination angles, and the 
error bars show the standard deviation. The dashed hues show the input value to the RT 
code. 
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Fig. 3. — Fraction of silicate mass in crystals {Left) and mass-averaged grain size (Right) for 
(a) Herbig Ae star, (6) T Tauri star and (c) Brown Dwarf. The spectra were fitted without 
carbon among the fitted dust species, however carbon is always present in the 2D RT models. 
The symbols show the averaged value over the different inclination angles, and the error bars 
show the standard deviation. The dashed hues show the input value to the RT code. 
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Fig. 4. — The effect of the fitted wavelength range on the results of the T Tauri spectra, 
fitted by the TLTD method (Left) and by the TT method (Right). The symbols show the 
averaged value over the different inclination angles, and the error bars show the standard 
deviation. The dashed lines show the input value to the RT code. 
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Fig. 5. — Normalized crystallinity as a function of mass- averaged grain size for the best 

configuration of each method. A wavelength range of 7-14 fim was used for the ST, TT 
and COS methods, while we used 7-17/Lim for the TLTD method. The spectra were fitted 
without carbon using the ST and TT methods while for the TLTD and COS methods carbon 
was included in the fit. 
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Fig. 6. — Normalized crystallinity as a function of mass- averaged grain size for all T Tauri 
spectra fitted without 6 //m size grains. A wavelength range of 8-13 /xm was used for all 
methods. Carbon was included in the fit only in the case of the TLTD method. 
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Fig. 7. — Mass-averaged grain size as a function of total crystallinity derived with the TLTD 
method. Both statistics are normahzed to the input value of the 2D RT code (indicated by 
the star in the middle of the figures). The fitted wavelength range was 8-13 /xm, while the 
applied signal to noise ratios were a), 10, h) 100, c) 1000. 



-40- 



1 














•■'?.>.:•"■:.■•.. ■ 




•■m. 


(o) 


• 


(b) 


(c) 





0.1 1.0 0.1 1.0 0.1 1.0 

<o>/<a > <a>/<a > <o>/<o > 

inp ' inp ' inp 



Fig. 8. — Mass-averaged grain size as a function of total crystallinity derived with the TLTD 
method. Both statistics are normahzed to the input value of the 2D RT code (indicated by 
the star in the middle of the figures). The fitted wavelength range was 7-17 /xm, while the 
applied signal to noise ratios were a) 10, b) 100, c) 1000. 
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Fig. 9. — (a): Contribution of different annuli to the total flux at 10 //m in the TLTD fit of 
a Herbig Ae disk model (Left) and a Brown Dwarf disk (Right) from the BS model series at 
an inclination angle of 45°. The grey area marks the region from which 70% of the total 
flux comes from, (b): Radial distribution of the temperature in the disk atmosphere. It has 
been calculated for the 2D RT model (filled circles) , the results of the TLTD code (hollow 
diamonds), the result of the TT method (dashed line) and the ST method (dot-dashed line), 
(c): Radial distribution of the temperature in the continuum. The symbols are the same as 
in the figure in the middle. 
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Fig. 10. — Spcctnim of a T Tauri model in the AS series at an inclination angle of 45° fitted 
by a) the ST method, b) the TT method and c) the TLTD method. The horizontal line 
shows the wavelength range of the fit. 




Fig. 11. — Contours of log{}(^) for the TT method during the fitting of the T Tauri model 
in the AS series at an inchnation angle of 45°. The black star marks the global minimum 
and the contour lines are separated by 0.5. 
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Fig. 12. — Effects of different parameters on the 10 /xm features, (a) Spectra of ffared disks 
around different types of stars at an inclination angle of 45°. (6) Spectra of T Tauri stars 
with different disk geometries (flared, moderately flared and flat) at an inclination angle of 
45°. Spectra of T Tauri stars with a flared disk but at different inclination angles. 
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Dust species 



1 : Am. Sil. (Olivine) 0.1 

2 : Am. Sil. (Olivine) 1 .5 /im 
5 : Am. Sil. (Olivine) 6.0 /j.m 

4 : Am. Sil. (Pyroxene) 0.1 /im 

5 : Am. Sil. (Pyroxene) 1.5 fj.m 

6 : Am. Sil. (Pyroxene) 6.0 /j.m 

7 : Forsterite 0.1 /im 

8 : Forsterite 1 .5 fim 

9 : Forsterite 6.0 fim 

10 : Enstotite 0.1 /u.m 
1 1 : Enstotite 1 .5 fj.m 
1 2 : Enstotite 6.0 

1 3 : Am. Silica 0. 1 //m 

14 : Am. Slllco 1.5 fj.m 

15 : Am. Silica 6.0 ^m 

16 : Am. Carbon 0.1 /.tm 

17 : Am. Carbon 1.5 fj.m 

18 : Am. Carbon 6.0 fim 



Fig. 13. — a (for the definition, see Sec. l4.4l) as a quality indicator of how well the optical 
data of a specific dust species can be reproduced by the rest. The curve of each dust 
species was fitted using all the other dust species in the wavelength range of (a) 8-13 /xm, 
(6) 7-14 /xm, (c) 7-1 7 /im and (d) 5-35 /xm . 



Table 1. Main parameters of the 2D RT disk models. 





Herbig Ae star 


T Taur star 


Brown Dwarf 




30 


1.4 


0.035 


T. [K] 


9500 


4000 


2500 


[Mq] 


3 


0.55 


0.08 


Rin [AU] 


0.4 


0.087 


0.014 


Rout [AU] 


400 


180 


87 



Table 2. Overview of dust species used. For each component we specify its lattice 
structure, chemical composition, shape and reference to the laboratory measurements of 
the optical constants. For the homogeneous spheres we used Mie theory to calculate the 
opacities. For the inhomoge neous spheres, we used the distribution of hollow spheres 
(DHS: lMin et al.ll2005al ). to simulate grain deviating from perfect symmetry. 



# 


Species 


State 


Chemical 
formula 


Shape 


Ref 


1 


Amorphous silicate 


A 


MgFeSi04 


Homogeneous sphere 


(1) 




(Olivine stoichiometry) 










2 


Amorphous silicate 


A 


MgFeSi206 


Homogeneous sphere 


(1) 




(Pyroxene stoichiometry) 










3 


Forsterite 


C 


Mg2Si04 


Hollow sphere 


(2) 


4 


Clino Enstatite 


C 


MgSiOs 


Hollow sphere 


(3) 


5 


Silica 


A 


SiOs 


Hollow sphere 


(4) 


6 


Amorphous Carbon 


A 


aC 


Homogeneous sphere 


(5) 



References. 



Dorschner et al. 


( 


1995 


1); (2 


Henning & Mutschkel ( 


1997) 



Servoin &; Piriou 


1973); 


(5) 


Preibisch et al. 


( 


1993) 



-47- 



Table 3. Result of the spectral decomposition of the spectrum of a Herbig Ae disk with an 
inchnation of 45°. Fitted wavelength domain : 7-14 //m. Similar tables for all fits can be 

found in the electronic edition of the Journal. 



Name 


Size 


M. Fraction 


M. fraction 


M. Fraction 


M. fraction 


M. fraction 




[^m] 


Real [%] 


cos [%] 


ST [%] 


TT [%] 


TLTD [%] 


wiivine yj\ j 


n 1 

U. 1 


lO.U 


u.u 


Q A 




lo.o 


wilvine J 




7 n 


-Ly.o 


n n 

U.U 


U.U 


O.C5 




o.u 


n 


U.U 


U.U 


U.U 


A "i 


Pyroxene (A) 


A 1 

U.i 


io.U 


U.U 


i.O 




1A A 


Pyroxene (A) 


1.5 


7.0 


14.8 


0.0 


0.0 


7.4 


Pyroxene (A) 


6.0 


2.0 


0.0 


0.0 


0.0 


2.8 


Forsterite (C) 


0.1 


2.0 


0.1 


1.7 


2.5 


1.9 


Forsterite (C) 


1.5 


1.0 


0.0 


0.0 


1.4 


1.0 


Forsterite (C) 


6.0 


0.0 


28.5 


0.0 


0.0 


0.1 


Enstatitc (C) 


0.1 


2.0 


0.0 


1.6 


2.6 


1.9 


Enstatite (C) 


1.5 


1.0 


1.7 


2.2 


0.1 


0.9 


Enstatite (C) 


6.0 


0.0 


11.2 


3.8 


7.6 


0.9 


Silica (A) 


0.1 


0.0 


0.0 


0.0 


0.0 


0.0 


Silica (A) 


1.5 


0.0 


0.0 


0.0 


0.0 


0.0 


Silica (A) 


6.0 


0.0 


10.0 


0.0 


0.0 


0.0 


Carbon (A) 


0.1 


46.0 


14.0 


1.6 


0.0 


43.8 


Carbon (A) 


1.5 


0.0 


0.0 


0.0 


0.0 


0.0 


Carbon (A) 


6.0 


0.0 


0.0 


78.2 


26.7 


0.0 


Average grain size 




1.00 


1.5 


0.10 


0.10 


1.35 


Total crystallinity 




6.00 


41.5 


9.30 


14.20 


6.79 


Tf [K] 








380.0 


190.0 




Tc [K] 








380.0 


615.0 




Trim(Rin) [K] 












1408.0 


Tatm(lAU) [K] 












755.6 


T^id(lAU) [K] 












171.9 


Prim 












-1.39 


Patm 












-0.35 


Pmid 












-0.10 


Reduced 






301.6 


2697.0 


33.19 


0.051 



